function energy2

%  comparison of energy using verlet and RK4
%  y'' = f(t,y)  with   y(0) = y0

% clear all previous variables and plots
clear *
clf

% parameters for calculation
n=8000;
t0=0; 
y0=[1 0];
tmin=3987;
tmax=4000;

%  calculate numerical solutions
t=linspace(t0,tmax,n);
h=t(2)-t(1);

%  RK4
[y_rk4,v_rk4]=rk4('de_f',t,y0,h,n);
for i=1:n
	e_rk4(i)=0.5*y_rk4(i)^2+0.5*v_rk4(i)^2;
end;

%  Verlet
[y_v,v_v]=verlet('de_f',t,y0,h,n);
for i=1:n
	e_v(i)=0.5*y_v(i)^2+0.5*v_v(i)^2;
end;

% get(gcf)
% set(gcf,'Position', [1041 771 548 229]);

% plot energy
subplot(1,2,1)

% plot solutions
te(1)=tmin;  te(2)=tmax;
h(1)=0.5; h(2)=0.5;
plot(te,h,'-k')
hold on
plot(t,e_v,'-r','LineWidth',1)
plot(t,e_rk4,'-b','LineWidth',1)
legend(' Exact',' Verlet',' RK4','Location','East');

% define axes used in plot
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('Energy','FontSize',14,'FontWeight','bold')
axis([tmin tmax 0.0 0.6001])

% have MATLAB use certain plot options (all are optional)
set(gca,'ytick',[0 0.1 0.2 0.3 0.4 0.5]);
set(gca,'xtick',[tmin tmax]);
axis([tmin tmax 0 0.52])
box on
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold');  
hold off

ic=1;
for it=tmin:tmax
	y4(ic)=y_rk4(it);  v4(ic)=v_rk4(it);
	yv(ic)=y_v(it);  vv(ic)=v_v(it);
ic=ic+1;
end;

tt=linspace(0,2*pi,100);
for id=1:100
	ye(id)=cos(tt(id));  ve(id)=-sin(tt(id));
end;

% plot phase plane
subplot(1,2,2)

% plot solutions
plot(ye,ve,'-k')
hold on
plot(y4,v4,'-b','LineWidth',1)
plot(yv,vv,'-r','LineWidth',1)

% define axes used in plot
axis equal
axis([-1.1 1.1 -1.1 1.1])
xlabel('y-axis','FontSize',14,'FontWeight','bold')
ylabel('v-axis','FontSize',14,'FontWeight','bold')

% have MATLAB use certain plot options (all are optional)
set(gca,'xtick',[-1 0 1]);
set(gca,'ytick',[-1 0 1]);
box on
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14);  
hold off



% right-hand side of DE
function z=de_f(t,y)
z=[y(2) -y(1)];


% RK4 method
function [ypoints, vpoints] =rk4(f,t,y0,h,n)
y=y0;
ypoints=y0(1);
vpoints=y0(2);
for i=2:n
	k1=h*feval(f,t(i-1),y);
	k2=h*feval(f,t(i-1)+0.5*h,y+0.5*k1);
	k3=h*feval(f,t(i-1)+0.5*h,y+0.5*k2);
	k4=h*feval(f,t(i),y+k3);
	yy=y+(k1+2*k2+2*k3+k4)/6;
	ypoints=[ypoints, yy(1)];
	vpoints=[vpoints, yy(2)];
	y=yy;
end;

% RK2 method
function [ypoints, vpoints] =rk2(f,t,y0,h,n)
y=y0;
ypoints=y0(1);
vpoints=y0(2);
for i=2:n
	k1=h*feval(f,t(i-1),y);
	k2=h*feval(f,t(i),y+k1);
	yy=y+0.5*(k1+k2);
	ypoints=[ypoints, yy(1)];
	vpoints=[vpoints, yy(2)];
	y=yy;
end;


% Verlet method
function [ypoints, vpoints] =verlet(f,t,y0,h,n)
ypoints=y0(1);
vpoints=y0(2);
ya=y0(1);
va=y0(2);
for i=2:n
	y=ya+h*va-0.5*h*h*ya;
	v=va+0.5*h*(-y-ya);
	ya=y;  va=v;
	ypoints=[ypoints, y];
	vpoints=[vpoints, v];
end;